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Abstract 



We define the notion of effective stiffness and sliow that it can used to build sparsifiers, 
algorithms that sparsify linear systems arising from finite-element discretizations of PDEs. In 
particular, we show that sampling O(nlogn) elements according to probabilities derived from 
, effective stiffnesses yields an high quality preconditioner that can be used to solve the linear 

' system in a small number of iterations. Effective stiffness generalizes the notion of effective 

, resistance, a key ingredient of recent progress in developing nearly linear symmetric diagonally 

dominant (SDD) linear solvers. Solving finite elements problems is of considerably more interest 
O . than the solution of SDD linear systems, since the finite element method is frequently used to 

numerically solve PDEs arising in scientific and engineering applications. Unlike SDD systems, 
which are relatively easy to precondition, there has been limited success in designing fast solvers 
^ ' for finite element systems, and previous algorithms usually target discretization of limited class 

, of PDEs like scalar elliptic or 2D trusses. Our sparsifier is general; it applies to a wide range 

£^ ■ of finite-element discretizations. A sparsifier does not constitute a complete linear solver. To 

construct a solver, one needs additional components (e.g., an efficient elimination or multilevel 
scheme for the sparsified system). Still, sparsifiers have been a critical tools in efficient SDD 
solvers, and we believe that our sparsifier will become a key ingredient in future finite-element 
solvers. 



1 Introduction 



^ I We explore the sparsification of finite element matrices using effective stijjness sampling. The 

goal of the sparsification is to reduce the number of elements in the matrix so that it can be 
easily factored and used as a preconditioner for an iterative linear solver. We show that sampling 
non-uniformly 0(n log n) elements produces a matrix that is with high probability spectrally close 
to the original matrix, and therefore an excellent preconditioner. The sampling probability of 
an element is given by the largest generalized eigenvalue of the element matrix and the effective 
stiffness matrix of the element. 

Effective stiffness generalizes the notion of effective resistance, a key ingredient in much of the 
recent progress in nearly optimal symmetric diagonally dominant (SDD) linear solvers [9l \2\ [10]. 
Solving finite elements problems is of considerably more interest than the solution of SDD linear 
systems, since the finite element method is frequently used to numerically solve PDEs arising in 
scientific and engineering applications. 

Unlike SDD systems, which are relatively easy to precondition, there has been limited success in 
designing fast solvers for finite element systems. Efforts to generalize combinatorial preconditioners 
to matrices that are not weighted Laplacians followed several paths, and started long before recent 
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progresses. Gremban showed how to transform a hnear system whose coefficient matrix is a 
signed Laplacian to a hnear system of twice the size whose matrix is a weighted Laplacian. The 
coefficient matrix is a 2-by-2 block matrix with diagonal blocks with the same sparsity pattern 
as the original matrix A and with identity off-diagonal blocks. A different approach is to extend 
Vaidya's construction to signed graphs [3|. The class of symmetric matrices with a symmetric 
factorization A = UU'^ where columns of U have at most 2 nonzeros contains not only signed 
graphs, but also gain graphs, which are not diagonally dominant [3]; it turns out that these 
matrices can be scaled to diagonal dominance, which allows graph preconditioners to be applied 
to them [7]. 

The matrices that arise in finite-element discretization of elliptic partial differential equations 
(PDEs) are positive semi-definite, but in general they are not diagonally dominant. However, 
when the PDE is scalar (e.g., describes a problem in electrostatics), the matrices can sometimes 
be approximated by diagonally dominant matrices. In this scheme, the coefficient matrix A is 
first approximated by a diagonally-dominant matrix D, and then G^) is used to construct the 
graph of the preconditioner B. For larg matrices of this class, the first step is expensive, but 
because finite-element matrices have a natural representation as a sum of very sparse matrices, the 
diagonally-dominant approximation can be constructed for each term in the sum separately. There 
are at least three ways to construct these approximations: during the finite-element discretization 
process [5], algebraically [1], and geometrically [18]. A slightly modified construction that can 
accommodate terms that do not have a close diagonally-dominant approximation works well in 
practice [l]. 

Another approach for constructing combinatorial preconditioners to finite element problems is 
to rely on a graph that describes the relations between neighboring elements. This graph is the 
dual of the finite-element mesh; elements in the mesh are the vertices of the graph. Once the 
graph is constructed, it can be sparsified much like subset preconditioners. This approach, which 
is applicable to vector problems like linear elasticity, was proposed in |14j : this paper also showed 
how to construct the dual graph algebraically and how to construct the finite-element problem 
that corresponds to the sparsified dual graph. The first effective preconditioner of this class was 
proposed in [6]. It is not yet known how to weigh the edges of the dual graph effectively, which 
limits the applicability of this method. However, in applications where there is no need to weigh 
the edges, the method is effective |15] . 

Our theory of effective stiffness sampling is an extension of the theory of effective resistance 
sampling. It is applicable to a wide range of finite element discretizations. But our sparsifier 
is not yet a complete algorithm for solving finite-element systems. We discuss the remaining 
challenges in Section El Nevertheless, we our results constitute a useful technique that should lead 
to fast finite-element solvers. A similar evolution gave rise to the fastest SDD solvers: Speilman 
and Srivastava's theory of effective resistance sampling |16] did not immediately lead to efficient 
algorithm, but the follow-up work of Koutis et al. turned it into very efficient algorithms [9l IIP]. 
The techniques used by the authors of [9l, [TO] to solve SDD systems do not trivially carry over to 
finite element matrices. For example, their constructions rely on low-stretch trees, a concept that 
does not have a natural extension for finite element matrices. But we expect such extensions to 
be developed in the future. 
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2 Preliminaries 



2.1 Sums of Random Matrices 



Approximating a matrix using random sampling can be viewed as a particular case of sums of ran- 
dom matrices. In the last few years there has been significant literature on showing concentration 
bounds on such sums [13\ [12\ [T7] . We use the following bound from [11] . 



Theorem 2.1. UH Theorem 1.1] Let Ai,A2, . 
that the AiS are real and symmetric with \\ E{A. 
and let M = il(7log(7/e^)/e^). // every value in the support of Ai has rank at most M then 



be i.i.d matrix-valued random variables. Assume 
2 < 1 and \\Ai\\2 < 7 almost surely. Let < e < 1 



Pr 



1 *^ 



1 



poly(M) 



2.2 Generalized eigenvalues, analysis of iterative methods and spar sificat ion 
bounds 

A well known property of many iterative linear solver, including the popular conjugate gradient 
and the theoretically convenient Chebyshev iteration, is that their convergence rate depends on 
the distribution of the eigenvalues of the coefficient matrix (its spectrum). The rate depends on 
how much the spectrum is clustered, but it is hard to form a concise bound. A simple and useful 
theoretical bound for symmetric positive semidefinite matrices depends only on the ratio between 
the largest and smallest eigenvalue. When using preconditioned methods convergence is governed 
by the generalized eigenvalues. 

Definition 2.2. Given two matrices A and S in M with the same null space N, a finite generalized 
eigenvalue A of {A, B) is a scalar satisfying Ax = \Bx for some x N. The generalized finite 
spectrum A{A, B) is the set of finite generalized eigenvalues of (^4, B), and the generalized condition 
number B) is 

maxA(^,^) 
mmA{A,B) ' 

We define the trace of {A, B) (denoted by Tr(yl, B)) as the sum of finite generalized eigenvalues of 
{A,B). 

(This definition can be generalized to the case of different null spaces, but this is irrelevant for 
this paper.) We will denote by A.{A) the set of finite non-zero eigenvalues of A (which is equal to 
A(^, Pa), where Pa is a projection onto the range of A). 

We are mainly interested in bounds on the smallest and largest generalized eigenvalues (which 
we denote Amin(', •) and Amax(') •) respectively), since they tell us two important properties on the 
pair (A, B). First, for every unit norm vector x we have 

Amm(^, B) ■ x'^Bx < x'^Ax < Amax(^, B) ■ x'^Bx . 

Second, when B is used as a preconditioner for A, a vector x satisfying — < e||^"'"&m is 

found in at most 0{y^ k{A, B) ■ log(l/e)) iterations where = x^Ax. 

In many cases it is easier to reason about non- generalized eigenvalues. The following result 
from [il relates generalized eigenvalues with regular eigenvalues of a different matrix. 
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Lemma 2.3. Let A = UU'^ and B = VV'^ , where U and V are real valued with the same number 
of rows. Assume that A and B are symmetric, positive semidefinite and {A) = (B). We have 

A{A,B) = {V+U) 

and 

A{A,B) = . 

In these expressions, S(-) is the set of nonzero singular values of the matrix within the parenthe- 
ses, Tj^ denotes the same singular values to the £th power, and denotes the Moore-Penrose 
pseudoinverse ofV. 

2.3 Effective resistance sampling 

Recent progress of in fast SDD solvers [HI EJ [10] is based on effective resistance sampling, first 
suggested in [16j. Solving SDD systems can be reduced to solving a Laplacian system. Given a 
weighted undirected graph G = {[n],E, w) the Laplacian L is given hy L = D — A where A is the 
weighted adjacency matrix Aij = Wij and D is the diagonal matrix of weighted degrees given by 
Da = Ylij^i'^ij- The effective resistance Re of an edge e = (u, is given by 

Lie — ipu ^v) L (c-u 6-^) 

where and Cy are identity vectors and L~^ is the Moore-Penrose pseudoinverse of L. The quantity 
is named effective resistance because Re is equal to the potential difference induced between u and 
V when a unit of current is injected at u and extracted at v, when G is viewed as an electrical 
network with conductances given by w. 

Spielman and Srivastava [16j showed that sampling sufficiently enough edges, where the prob- 
ability of sampling an edge is proportional to WeRe yields an high-quality sparsifier for G, which 
can be translated to an high-quality preconditioner. Koutis et al. [9l El |T0] show that even crude 
approximations to the accurate effective resistances suffice, and they show how such an approxi- 
mation can be computed efficiently. The asymptotically fastest solver |10j solves an n-by-n SDD 
linear system in time 0(m log n log (1/e)) where m is the number of non-zeros in the matrix and e 
is the accuracy of the solution. 

3 Finite Element Matrices and Their Factored Form 

A finite element discretization of a PDE usually leads to an algebraic system of equations Kx = b. 
The matrix K has certain properties that stem from the PDE and the specifics of how it was 
discretized. To make our results more general and easier to understand by a wide audience, we use 
the algebraic-combinatorial formulation developed in [14] rather than a PDE-derived formulation. 

The matrix K £ M"^"- is called a stiffness matrix, and it is a sum of element matrices, K = 
X^^i ^e- Each element matrix Ke corresponds to a subset of the domain called a finite element. 
The elements are disjoint except perhaps for their boundaries and their union is the domain. We 
assume that each element matrix Ke is symmetric, positive semidefinite, and zero outside a small 
set of Ue rows and columns. In most cases Ue is uniformly bounded by a small integer. We denote 
the set of nonzero rows and columns of Ke by Me. We denote the restriction of a matrix A to 
indices / by A{L), and denote the Ke = Ke{N'e). Ke is the essential element matrix of e. Typically, 
in finite element discretizations both the stiffness matrix [K) and the essential element matrices 
{Kes) are singular. We denote the dimension of the null space of Ke by de = dim(null(-fCe)) and 
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the rank of by = Ue — dg. For simplicity, we will assume the rank (and dimension of null 
space) of all the elements is the same and equal to r (d for null space dimension). The results can 
be easily extended for non uniform element ranks. The null space of K is denoted by N and we 
assume that its dimension is d as well. 

Most, if not all, theoretical results on sampling matrices are on rectangular matrices, and their 
usefulness rely on the aspect ratio being high. Finite element matrices, according to our definition, 
are square. Luckily, K can be written as K = F'^F where F has more rows than columns. The 
key is obtaining a factored form of = FJ Ff> where is r x n. We can then write 



and K = F'^F. Many finite-element discretization techniques actually generate the element ma- 
trices in a factored form. Even if the elements are not generated in a factored form, a factored 
form can be easily computed. One way to do so is using the eigendecomposition = Vf^eVj ■ 
Define F^. = E^^Vg^ where V^. is obtained by taking the r columns of 14 associated with non-zero 
eigenvalues, and let Ff^ be obtained by expanding the number of columns of F^. to n by adding zero 
columns for columns not in Me- It is easy to verify that K^, = FJ Ff> and that is r x n. 

Typically, the element matrices are compatible with the null space of K meaning that the 
null space of is the restriction of the null space of K to Me- When all the element matrices 
are compatible with the null space of K, the matrices involved have useful properties p[4] that we 
use in theorems as needed. In the Appendix we elaborate on the issue of null-space compatibility 
and explain why finite element matrices typically have the properties assumed in the theoretical 
analysis. 

The main property we assume is that the rank deficiency of the factor F is minimal. 

Definition 3.1. A matrix F G ]^»"X" Y^g^ minimal rank deficiency if every set of n — dim(null(-F)) 
columns of F is independent. 

Note that if the rank deficiency of F is minimal then every leading / x / minor of K is non- 
singular, as long as I < n — d. 

4 Effective Stiffness of an Element 

We now define the effective stiffness of an element. The stiffness matrix of an element describes 
the physical properties (elasticity, electrical conductivity, thermal conductivity, etc) of a piece of 
material called an element by showing how that piece of material responds to a load (current, 
mechanical force, etc) placed on the element. The effective stiffness matrix shows how the entire 
structure responds to a load that is placed on one element. Intuitively, if the stiffness matrix and 
the effective stiffness matrix of an element are similar, the element is important; removing it from 
the structure may significantly change the behavior of the overall structure. On the other hand, if 
the effective stiffness element has a much larger norm than the element matrix, then the element 
does not contribute much to the strength (or conductivity) of the overall structure, so it can be 
removed without changing much the overall behavior. 

Algebraically, the effective stiffness matrix of e is obtained by eliminating from K all columns 
not associated with e. 
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Definition 4.1. Let K be obtained from K by an arbitrary symmetric reordering of the row 
and columns of K such that the last Ue rows and columns of K are J\fe and they are ordered in 
ascending order (i.e., the ordering in K of the columns in J\fe is consistent with their order in K). 
Suppose that K is partitioned 

^-{Rf, K22 ) 

where Ki G M("-"e)x(n-ne)^ G r("-"-)x"- and K22 G M"-^"<=. If Ku is non singular we say 
that element e is supported . The effective stiffness Se of element e is 

^ _ \ K22 - Kj2K^i K12 e is supported 
I Kf, otherwise 

It is easy to verify that the the effective stiffness is well defined in the sense that any ordering 
that respects the conditions of the definition gives the same Sg- Note that if the factor F has 
minimal rank deficiency then all elements are supported. 

Before proceeding to discuss effective stiffness sampling, and stating our main result, we first 
show that indeed effective stiffness generalizes effective resistance by showing that effective resis- 
tance is a particular case of effective stiffness. 

The Laplacian of a weighted graph G = {[n],E,w) is, in fact, a finite element matrix per our 
definition in section [3l Given an edge e = {u,v) define Kg = We{eu — ev){eu — e^)^ . It is easy to 
verify that L = Yle<=E ^e- well-known that if the if the graph is connected the rank deficiency 
is d = 1 and the null space N is all-ones vector. L can also be written in factor form L = F"^ F 
where F G M'^I^I^L Each edge e = {u,v) correspond to row in F given by Fg = ^/wl{eu — e^)'^ . 
If the graph is connected the factor F has minimal rank deficiency, so all elements (edges) are 
supported. 

Lemma 4.2. Let F be the factor of a graph G. If G is connected then F has minimal rank 
deficiency. 

Proof. Suppose there is a non-independent size n — 1 subset of F's columns. Let F be a reorder 
of F's columns such that those n — 1 columns are the first n — 1 columns. Since the first n — 1 
columns of F are linearly dependent there is a vector x 7^ such that 

The vector ( is not spanned by 1„. This contradicts our assumption that G is connected 
since the null space of the Laplacian of a connected graph is always span{l„}. □ 

Simple calculation shows that 5612 = and (ei — e2)'^Se{ei — 62) = Re^- This implies that 
Se — Rg {c-u 6?))(c-!i ey) • 

Graph sparsification by effective resistance [16] and near- linear time linear solvers [9l[2l[T0] relay 
on sampling edges with probability relative to WeRe- It is easy to verify that WeRe = ^ina.x{Ke, Se)- 
We call the quantity Amax(-^e, •S'e) the leverage of element e. Our main result shows that sampling 
probabilities should be relative to the leverages for general finite element matrices, and not only 
for Laplacians. 
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5 Effective Stiffness Sampling 

This section defines the leverage of an of element and shows that sampling probabilities that are 
relative to the element leverages is a good choice. 

Definition 5.1. The leverage of an element e is 
The total leverage of K is 

m 
e=l 

Note 5.2. The term leverage arises from the connection between effective resistance and statistical 
leverage that was noted by Drineas and Mahoney in [8]. 

Lemma 5.3. Let K = F^F = ^^j^ K^, be an n-by-n finite element matrix. We have {n — d)/r < 
tk < n — d. 

We prove Lemma 15.31 later this section and first state the main theorem of this paper. The 
main theorem shows how to use the leverages to sample finite element matrices. 

Theorem 5.4. Let K = F^F = Xl^Li Ke be an n-by-n finite element matrix. Assume that the 
factor F has minimal rank deficiency and null(S'e) = null(-K'e) for every element e. Let 



and let Ti, . . . , Tm be a i.i.d random matrices defined by 

where Ji, . . . , Jm ore random integers between 1 and m which takes value e with probability p^. In 
other words, Ti is a scaled version of one of the K^s, selected at random, with a scaling that is 
proportional to the inverse of p^. For M = ^{TK^og{TK)) (and M = Q.{n\og{n)) ) , we have 

/ 1 ^ \ 1 

To prove Lemma 15.31 and Theorem 15.41 we need the following Lemma. 

Lemma 5.5. Let U G i^^f xn ^q^ij-^x whose columns form an orthonormal basis of range{F) . 

Let Ue € M'^^" be the rows ofU corresponding to element e. Assume that the factor F has minimal 
rank deficiency and null(Se) = null(-fCe)- I'he set of non-zero eigenvalues (including multiplicity) 
of UeU^ and the set of finite generalized eigenvalues of {Ke,Se) are the same. In particular, 

and 

Tr{UeUl) = Tr{k„Se). 



7 



Proof. We first show that we can prove the lemma by showing that it holds for a particular U. An 
arbitrary orthonormal basis V is related to ?7 by y = UZ, where Z is an n-by-n unitary matrix. 
In particular, Ve = UeZ {Ve are the rows of V corresponding to element e) so VeVj' = UeZZ^Uj = 
UeU^ . We obtain U from the QR factorization of F = UR and set U to be the first n — d columns 
of U, where F is obtained from F by reordering the columns in Me to the end (consistently with 
their ordering in F). 

The last rig columns of F are J\fe, and Fg is non-zero outside the indices of AAe.This implies that 

Fe = [ Orx(n—ne) -^e ] 
Ue = [ Orx(n—ne) ] 

where Ue, Fe £ W^"^" . Let us write 



R 



Rii Ri2 
R22 



where Rn G m("-"'^)x("-"-), R^o e m("-"'=)x"'= and R22 e ^"■^x"^ Let us write K = F'^ F and 



where Ku E M("-"e)x{n-ne)^ j;^^^ G M^"""-^)^"- and K22 e M"-^"-. Since i? is the i?-factor of F 
and K = F'^F it is also the Cholesky factor of K. It also implies that R22R22 is equal to the Schur 
complement 

-^22-^22 = -^22 — Ki2^ll ^12 = Sg ■ 

The minimal rank deficiency of F implies that that the bottom d rows of R and i?22 are zero. 
Let ^22 e m('^<=-'^)x"<= be the first He - d rows of i?22- It is still the case that -^^"2-^22 = -^e. We 
have F = UR, so F^ = UeR22 = UeR22 which implies that Fg = UeR22- Applying Lemma 12.31 we 
find that 

A{Ke,Se) = A(FjFe,i?i2^22) 

= J:'{{R^2rF. 

= ( {R22)^R22Ue 



The minimal rank deficiency of F implies that -R22 is full rank, so R^2 is a full rank matrix with 
more rows than columns (or equal), so (-^22)^-^22 — ^rie- This implies that {R22)~^^22^T — so 



A(^e,5e) =S2(C7J) 



T,'^{Uj') is exactly the set of non-zero eigenvalues of UeU'^ ■ Therefore, the non-zero eigenvalues of 
UelJj are exactly the finite generalized eigenvalues of {Ke,Se), so 

■^max ( t^e f^e ) — '^max(-^e) 'S'e) — 

and 

Tr{UeUj)=Tl{ke,Se). 

□ 
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We can now prove Lemma 15.51 and Theorem 15.41 



Proof. ( of Lemma \5, 



m 

TK = Y.'^e 

e=l 

m 



— ^ ^ Amax(-^e; 'S'e) 
e=l 

m 

e=l 

m 

e=l 

m 

= Y.TT{U^Ue) 

e=l 

m 
i=l 

= Tr(C/^[/) = n-d. 

For each element the pencil {Ke,Se) has exactly r determined eigenvalues, so Amax(-^e) S'e) > 
Tr(Ere, S'e)/r. The lower bound follows. □ 

Proof. ( of Theorem \5.4\ ) We express the matrix Z^^lf i as a normal form 

M 



-^5:T, = (5Ff(5F) 



M 

where 5 G -^Mrxmr -g random sampling matrix and F is the factor of the stiffness matrix 
K = F^F. If we take 5 to be a block matrix with r x r blocks, its blocks defined by 

JjPe ^ Irxr if Ti = ^Ke 

[Orxr otherwise, 

then it is easy to verify that S indeed satisfies the identity above. Let F = UR be a reduced QR 
factorization of F. The minimal rank deficiency of F implies that that the bottom d rows of R are 
zero. Let R G be the first n — d rows of R, and U G i^mrx{n-d) ^j^g ^^.g^ n — d columns 

of U. It is easy to verify that F = UR and F^F = R^R. R^ is full rank, so {R^)+R^ = In. 
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Assume for now that null(5-F) = null(F). Applying lemma [2T3l we have 



1 ^' 



M 

i=l 

= KiR^R,{SFf{SF)) 

= -'{(suf) 

= k\SU) 

= KiiSUfiSU)). 

To bound k{{SU)^{SU)) with high probability we first bound \\{SU)'^ (SU) — Inxn\\2 with 
high probability. Define the i.i.d random matrices Yi, . . . , Ym by 

y, = pjIuIuj^ 

where Ug is the rows corresponding to element e in U. It is easy to verify that 

1 



M 

i=l 



The expectation of the l^'s is the identity matrix, 

M 



E(li) = J2^v{Ti=pJ^Kj)pfufUj 

M 
M 



U U = In-dxn-d 



and their 2-norm is bounded by 



|Fj||2 < rcL&-KPj^\\uJUj^^ 

= maxpj'^X^i,^{UjUj) 



maxn. ^Xma.AKj,Sj) 

maxp~^Tj 
j ^ 

-1 \ 



max 
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We showed that || E(^i)ll2 = ||-^nxn|| = 1 and that ||li||2 <TK<n, so we can apply Theorem 12. II to 
the YiS with e = 1/3 and 7 = tk- The application of the theorem shows that for M = ^{tk log(Ti^)) 
(and M = il(n log(n))), 



Pr 



1 ^ 
-VYi 



1=1 



-I < 1 

^ 3 ) - poly(M) 



This bounds ||(5C/) (5f/) — /nxnUg with high probability. Finally, we note that for every sym- 
metric metric A, if \\A - I\\2 <t <l then k{A) < Applying this to {SU)'^{SU) with t = 1/3 
we find that whenever \\{SU)'^{SU) — Inxn\\2 — 1/3 (which happens with high probability), we 
have k{{SU)^{SU)) < 2. 

We stiU have to prove that nuU(cSF) = nun(F). U is full rank so null(F) = null(?7i?) = null(i?). 
On the otherhand SF = SUR, so as long as SU is full rank we have null(5-F) = null(i?). SU is 
fuh rank if and only if {SUf{SU) is non-smg ular. When \\{SU)'^{SU) - Inxn\\2 < 1/3 it is also 
the case that (SU)'^ (SU) is non-singular. □ 



6 Sampling Using Inexact Leverages or Upper Bounds 

Theorem 15.41 shows that the sampling probabilities that are proportional to Tg are effective for 
randomly selecting a good subset of elements to serve as a preconditioner. In practice it may be 
possible to obtain only estimates for the true maximum eigenvalues. The following two general- 
izations of Theorem 15.41 show that even crude approximations or upper bounds of the leverages 
suffice, provided that the number of samples is enlarged accordingly. 

Theorem 6.1. For every element e let be (1 -|- 6) -approximations to Tg, that is 

\fe - Tel <5-Te 

We make the same assumptions and use the same notation as in Theorem \5.4\ except that the 
probabilities pe are now given by 

- 

Pe ~ • 

For M = Q{TK/3log{TKf3)), where (3 = j^, we have 

/ 1 \ 1 

Pr k(K, — V Ti) > 2 < — - . 

Proof. The proof is identical to the proof of Theorem 15.41 except that the bound on ||li||2 needs 
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to be modified as follows: 

||yi||2 < maxpJ^\\UjUj\\^ 



max 

j 



max p. Xma.AKj,Sj] 



max 



V \Z^i=l '« 



< max 



(1 - 5)Te 



m 



-1 



e=l 



□ 



Theorem 6.2. For every element e let fg he and upper hound on r^, and let fx = YlT=i^e- 
We make the same assumptions and use the same notation as in Theorem \5.4\ except that the 
prohabilities Pe are now given hy 

Te 

Pe = —■ 

TK 

For M = il{fK log(fi^)) we have 



ly(M) 

Proof. The proof is identical to the proof of Theorem 15.41 except that the bound on ||li||2 needs 
to be modified as follows: 

\\Yi\\^ < maxpf\\uJUj\\^ 

= m.a-Kp~^\raa.yi{UjUj) 
= max ^ Amax {Kj ,Sj) 

j 

(ffe 

= max I — 

i \\TK 
Tj 

< tk • max — 

j Tj 

< tk . 

□ 



7 A Condition-number Formula for The Leverages 

In this section we show that the leverage Te can also be defined in terms of the condition number 
of {K, K — Ke). This condition number is the one related to preconditioning K by removing only 
element e. 
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Theorem 7.1. Let K = F^F = X^e^i -^e be an n-by-n finite element matrix. Assume that the 
factor F has minimal rank deficiency. For every element e, if e is supported and null(i^ — Ke) = 
null(X) then 

k{K, K -Ke)-l 



k{K,K-K,) 



otherwise Te 



1. 



Ke and Te = Amax(i^e, Ke) 



1 and the claim holds. Assume 



Proof. If e is not supported then Se 
e is supported. 

We now argue that if rank(Er — Ke) < rank(i^) then Te = 1. Let K be obtained from K — Kg 
by an arbitrary symmetric reordering of the row and columns of K such that the last rig rows 
and columns of K are Me and they are ordered in ascending order (i.e., the ordering in K of the 
columns in Me is consistent with their order in K). Suppose that K is partitioned 



K 



Ku ^12 
-^12 ^22 



where Ki G ]R(n-"e)x(n-n,)^ ^ 
non-singular. This implies that iank{K — K, 



rank(^) = n — Ue 



K22 — ^12^11 ^12 is the Schur complement. It is easy to see that K22 



. Since e is supported Kn is 
rank(A'22 —-^12-^11^-^12) since 

Kj^K^^Ki2 = Se-ke. 

On the other hand, using similar observations we find that rank(i^) = n — Ue + rank(S'e). Prom 
rank(ivr — Ke) < rank(i^) we find that rank(S'e — Ke) < rank(S'e). Therefore there exists a vector 
X such that SeX 7^ but {Se — Ke)x = 0. That x is an eigenvector of {Ke,Se) corresponding to 
the eigenvalue 1 since we have KeX = SeX but SeX 7^ 0. All eigenvalues of {Ke,Se) are bounded 
by 1 so we found that AmaxC-f^^e, •S'e) = 1- 

We now analyze the spectrum of {K,K — Ke). Without loss of generality assume e = m. Let 
S G ]R(ni-i)rxmr defined as 



'5 [ -^(m—l)rxmr 0(m— l)rxr ] • 

It is easy to verify that K - Ke = {SF)^{SF). 

Let F = UR be a reduced QR factorization of F. The minimal rank deficiency of F implies 
that that the bottom d rows of R are zero. Let R G be the first n — d rows of R, and 

U e M"»^x("-'^) be the first n-d columns of U. It is easy to verify that F = UR and F'^F = R^ R. 
The matrix R'- is full rank, so {R^)^R^ = In. U has orthonormal rows so U^U = I {n-d) x {n-d)- 
Applying lemma [2^ we have 



K{K,K-Ke 



A{F^F,{SF)^{SF)) 
A{R^R, 

S2 {iR^)+F^S^) 
T? {iR^)+R^U^S^) 
T?{U^S^) 
A{U^S^SU) . 



Let T G 



be defined as 



T 







(m— 1) X (m— l)r 
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It is easy to verify that S^S = Imrxmr — T . We now have 

AiU^S^SU) = A{U^ {Imrxmr -T)U) 

= A{U^U -U^TU) 

= J^{I{n-d)x(n~d) - W^TU) . 

Let Ue be the bottom r rows of U. It is easy to verify that U^TU = U^Ue, so A{U^S^SU) = 
A{I — UjUe)- Let (A,x) be an eigenpair of Uj^Ue, that is UjUeX = Ax. We have 

(/ - UjUe)x = X- UjUeX = X-\X={1-\)X, 

so (1 — A, x) is an eigenpair of / — UjUe- UjUe is an order n — d matrix of rank r < n — (i so it is 
singular. UjUe is also positive semidefinite so all its eigenvalues are non-negative. The last three 
facts imply that Xma.x{I—Uj^Ue) = 1. On the other hand, clearly Amin(-^— t^Jf^e) = 1 — Amax(f^Jf^e)• 
Combining these two together we find that 

.(A-,A--A,) = .(/-C/JC/.) = ^3ld^. 

This implies that 

K{K,K-Ke) -^--^^^^-^ 

The non-zero eigenvalues oiUeUj are exactly the non-zero eigenvalues of UjUe, so Xma.x{UeUj) = 
Xuia,yi{UjUe)- U '\s & matrix whose columns form an orthonormal basis of range(-F), so according 
to Lemma 15.51 we have Xmax{UeUj') = Tg, which concludes the proof. □ 



8 Discussion and Conclusions 

The results in this paper do not constitute practical solver. What are the remaining challenges 
that need to be addressed to construct a complete solver? 

• Computing the leverages. Neither of the two formulas for the leverages can be computed 
more efficiently than solving the linear system itself. Theorems 16.11 and 16.21 show that an 
approximation or an upper bound of the true leverages suffice. Using the effective-stiffness 
formulas on an element e and the elements within some distance from it (rather than on the 
entire finite-element mesh) yields an upper bound fg > Tg. We believe that this upper bound 
would be effective for sparsifying large meshes. 

• Number of elements in the sparsified system. Theoretically, the number of elements in 
an n-by-n finite element matrix can be as large as Q{n'^), in which case 0(n log n) elements 
is a big improvement. In practice, there are typically only 0{n) elements, so sampling 
0(n log n) element is not an improvement. It is worth noting that elements are sampled 
with repetition so in practice fewer than 0(n log n) distinct elements are sampled. If the 
probabilities (leverages) are highly skewed then the number of element can sampled can even 
approach n. An illustrative, but rather useless, example is the following. Consider a finite 
element matrix with exactly n elements with leverage 1 and all other elements with leverage 
0. All samples will be inside the group of elements with leverage 1, so there will only be n 
distinct elements in the sample. The authors of [9l [10] used highly skewed probabilities to 
handle SDD matrices with only 0(n) non-zeros. 
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• Non-zeros in factor. Once the finite element matrix has been sparsified, the sparsified 
matrix has to be factored, or it can serve as a foundation for a multilevel scheme. The cost 
of factoring sparsified the matrix and the cost of each iteration of PCG depend mainly on 
the number of non-zeros in the factor (the fill-in) and not on the number of non-zeros in the 
sparsified matrix. To build an effective preconditioner using sampling, the sampling must 
be guided so that the sampled matrix will have low fill. For the sparsifier to be useful in a 
multilevel scheme, the sparsified matrix must be easy to coarsen (eliminate vertices, faces, 
or elements to obtain a small mesh on which the process can be repeated). 

The same issues prevented the initial theoretical results of [16] from immediately producing a 
fast algorithm. But a few years later fast algorithms based of effective resistance sampling were 
suggested. The extension of effective resistance to effective stiffness is not trivial. We should not 
expect the other techniques used in SDD solvers to trivially extend to finite-element matrices. 
For example, the first step in SDD solvers is forming a low-stretch tree. There is currently no 
equivalent combinatorial object for finite-element matrices. 

Hopefully, our first step will be followed by additional ones that will enable the construction 
of general and efficient finite-element solvers. 
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9 Appendix: Null Space Compatibility and Rigidity of Finite El- 
ement Matrices 

In this section we review facts from the algebraic-combinatorial theory of rigidity of finite-element 
matrices developed in |14j . and show that it implies typical finite-element matrices have the prop- 
erties assumed in the results. This section is not necessary for the validity of the theoretical results, 
but it does relate them to the actual application. 

Typically, the element matrices are compatible with the null space of K. We now define what 
this means exactly. 

Definition 9.1. Let A be an m-by-n matrix, let Za be the set of its zero columns. We define the 
essential null space of A (enull(^)) by 

enull(^) = {x : Ax = and Xj = for z G Za} ■ 
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Definition 9.2. Let N C be a linear space. A matrix A is called 'H- compatible (or compatible 
with N) if every vector in enull(A) has a unique extension into a vector in N, and if the restriction 
of every vector in N to Ma (setting indices outside Ma to zero) is always in enull(A). 

A particular discretization of a PDE yields element matrices {K^s) that are compatible with 
some well-known null space N, which depends on the PDE; a translation in electrostatics, transla- 
tions and rotations in elasticity, and so on. Furthermore, it is usually desirable that the stiffness 
matrix K be rigid with respect to N, which is equivalent to saying that the null space of K is 
exactly N. For example, for matrix of a resistive network (the Laplacian) elements are compatible 
with the span of the all-ones vector. The null space of the Laplacian is exactly the span of the 
all-ones (i.e., the matrix is rigid) if and only if the graph is connected. 

Lack of rigidity often implies that the PDE has not been discretized correctly, and it does 
not make sense to solve the linear equations. This is an important scenario to detect (see [15j), 
but it is not the subject of this paper. We will assume the matrix K is rigid with respect to the 
prescribed and well know null space. The prescribed null space typically (that is, for real-life finite 
element matrices) implies minimal rank deficiency, which has to be proved for each case. A simple 
technique (which the proof of Lemma 14.21 implicitlv uses) is based on the following lemma. 



// no 



Lemma 9.3. Suppose that K = F-^F G /f"X" Jig^g null space range(A^) where N G R"^°'. 
d X d suhmatrix of N is singular then F has minimal rank deficiency. 

Proof. First notice that null(-F) = null(i^) since null(F^) = range(F)"'". Suppose there is a set of 
n — d columns of F which are not independent. Let F be a reordering of the columns of F such 
that those n — d columns are first. There is a vector x G such that 



X 



0. 



Let be a reordering of the rows of A^ consistently with the reordering of the columns of F in 



F. The vector ( x" 

Ny=(- 



) is in the null space of F so there must exist a vector y ^ such that 

) ^ . This implies that the bottom d rows of N form a singular matrix. These rows 
are also rows of A^, which implies that A^ has a d x d singular submatrix, which contradicts our 
assumption. □ 

We already saw that this technique was used to show that the factor of a connected Laplacian 
has minimal rank deficiency. We now show another example: elastic struts in two dimensions. 
In [T3] it is shown that given a collection P = {pi}"^^of points in the plane, the null space of the 
rigid finite element matrix representing a collection of elastic struts between the points is spanned 
by the range of 

/I -yi\ 

1 Xl 

1 -y2 

1 X2 



N 



1 

Vo 



-yn 



The matrix A^ does not have singular 3-by-3 submatrix unless the points have some special proper- 
ties (like three points with the same x coordinate), which they typically do not have. Even if such 
a property is present, a slight rotation of the point set, an operation that does not fundamentally 
change the physical problem, will remove it. 
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The minimal rank deficiency of F implies that all elements are supported. We now show that 
null(5e) = null(-ftre) if K is rigid. To do so we need an additional definition and two lemmas 
from [14j . 

Definition 9.4. A matrix A is rigid with respect to another matrix B if for every vector in 
X G enul^A) there is a unique vector y G enull(i?) such that Xi = yi for all i G Na ^■N'b- The two 
matrices are called mutually rigid if they are rigid with respect to each other. 

Lemma 9.5. (Lemma 5.5 from Let N be a linear space, and let B be some matrix with no 
zero columns whose null space is N. Another matrix A is N-compatible if and only if A and B are 
mutually rigid. 

Lemma 9.6. (part of Lemma 3.7 from [I4]) Let A and B be n-by-n matrices of the form 



" ■ 


, B = 


' Bu 


B12 


A22 


B21 


B22 



Assuming that B has no zero columns, Bu is non-singular and A are mutually rigid we have 
nun(A22) = nun(S22 - B2iB{^^Bi2). 

It is easy to see that combining the last two Lemmas with minimal rank deficiency of F ensures 
that null(5'e) = null(/Ce) for every element e. 

Lemma 9.7. Let K = F^ F = ^^1^=1 n-by-n finite element matrix with null space 

N. Assuming that F has minimal rank deficiency, and that all element are N— compatible, then 
null(S'e) = null(Ke) for every element e. 
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